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THEOREnCAL  PREDICTION 
OF  VIBRATIONAL  CIRCULAR  DICHROISM  SPECTRA 
OF  R-GLYCERALDEHYDE,  R-ERYTHROSE,  AND  R-THREOSE 

n.  DEVELOPMENT  OF  A  PROCEDURE  TO  SCALE 
THE  FORCE  CONSTANT  MATRIX  EXPRESSED  IN  TERMS 
OF  INTERNAL  COORDINATES 


n.l  SCALING  PROCEDURES 


In  Pan  n  of  this  report  the  procedure  that  was  used  in  Part  I  to  scale  the  force  constant  matrix 
is  developed.  More  specifically  the  scaling  methods  used  are  discussed  and  the  programs  to 
implement  the  scaling  are  described  and  given.  The  harmonic  force  constant  matrix  gives  the 
second  derivative  of  the  energy  with  respect  to  the  coordinates  of  the  molecule.  The  force 
constant  can  be  expressed  in  two  common  ways.  First,  the  force  constant  matrix  can  be 
expressed  in  terms  of  Cartesian  coordinates  as 


Kij=d"Eldqidqj 


(1) 


where  is  a  Cartesian  coordinate  q,  =  x,,  q^  =  y,.  q,  =  z, . q.  ,  =  x  .  q.  .  =  y  .  and  q.  =  z  . 

where  n  is  the  number  of  atoms  in  the  molecule.  The  matrix  K  is  3n  x  3n.  Second,  the  force 
constant  matrix  can  be  expressed  in  terms  of  3n-6  internal  coordinates,  R;,  expressed  as  bond 
stretches,  bond  angle  bends,  dihedral  angle  torsion  or  other  modes  of  'motion.  This  force 
constant  matrix  is 


Fjj=d'EldRiBRj  (2) 

The  matrix  F  is  {3n-6)  x  (3n-6).  There  is  a  transformation  between  the  internal  and  Cartesian 
coordinates  given  by 


R=Bq  (3) 

where  R  and  q  are  column  vectors  whose  components  are  the  internal  and  Cartesian  coordinates. 
.Note  that  B  is  (3n-6)  x  (3n).  An  existing  computer  program  [1]  to  determine  B  given  R  and  q 
was  modified  and  named  bmat.f.  The  program  is  given  in  the  program  section  n.2;  in  addidon, 
in  Table  1  a  sample  setup  of  a  datafile  for  R-glycer^dehyde,  for  use  with  bmatf,  is  shown. 

The  following  relation  [2,3],  originally  shown  by  Pulay,  between  K  and  F  holds 

F=B-‘  KB-‘  -£<|>,B-‘  CB-‘  (4) 


where  4>j  is  the  column  vector  of  the  forces  expressed  in  internal  coordinates,  C'  is  the  second- 
order  transfmmatation  matrix  relating  the  Cartesian  and  internal  coordinates,  B  ‘  is  the  transpose 
of  B^  and  B*^  is  given  by 

B“*=(BinB'^r*Bm  (5) 

where  m  is  any  matrix  for  which  (B  m  B  '^)  is  not  singular.  [In  the  examples  studied  in  this 
report,  nonsingular  matrices  are  obtained  if  m  is  taken  to  be  the  identity  matrix.]  We  have 
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considoed  the  fosce  constant  matrices  at  an  qMifnued  geometry;  under  that  condition  equation  4 
becomes 

(6) 

A  program  named  matnnilLf  was  written  to  cany  out  Ac  matrix  multiplication  given  by  die 
pRwious  equations.  As  a  check  on  the  progntn  both  Gaussian  90  and  Gaussian  92  calculations 
on  optimu^  geometries  were  run  with  option  FREQ.  Hus  procedure  will  generate  bodi  force 
constant  matrices  K  and  F.  The  results  of  using  die  above  matrix  multiplication,  fw  the 
exanqiles  considered,  all  a{^  exaedy  with  the  results  obtained  frimi  the  Gaussian  calculations. 
The  program  matmultf  is  given  in  section  IL2. 

Next,  the  FORTRAN  program  matmultf  was  modified  to  allow  for  scaling  of  the  force 
constant  matrix,  F.  This  new  program  is  called  matmult2.f.  The  scaling  constants  are  input 
into  the  program  by  editing  matmult2.f.  The  resulting  scaled  F  matrix  is  converted  to  a  scaled  K 
matrix  which  is  then  used  as  input  to  the  CADPAC  program  to  carry  out  a  VQ>  calculation  of 
allowed  frequencies  of  vibrations  and  corre^nding  rotational  strengdis. 

IL2  SCALING  PROGRAMS 

In  this  section  3  FORTRAN  programs  are  reported. 

n.2.1  Program  Imiatf 

The  first  program  bmatf  determines  the  B  matrix  in  the  transformation  between  the  internal 
R  and  Cartesian  coordinates  q  where  R  is  a  (3n-6)-column  vector  and  q  is  a  (3n)-column  vector. 
The  data  file  used  must  be  set  for  each  molecule  that  is  considered.  Table  1  lists  a  sample  data 
file  for  R'glyceraldehyde. 

On  the  following  pages  a  listing  of  the  FORTRAN  program  bmatf  is  given. 


bnat.f 


.nf 

C»342  GEN  VZB  ANAL  PGM  USING  miSCM  GP  MATRIX  METHOD 
C 

C  THIS  IS  PROGRAM  NUMBER  1  OP  THE  COMPLETE  VIBRATIONAL  PACKAGE. 

C 

C  niAT  . . .  WIL.»(»}  B  MATRIX  ELEMENTS  POR  INTERNAL  COORDINATES 
C  iVERSIWO  JUL  28.  1977) 

C 

C  AUTHORSO  MIKE  PETERS(»1  AND  DOUG  MCINTOSH,  U  OP  T  CHBI  DEPT,  TORONTO 
C 

C  INPUTO 
C 

C  2  TITLE  CAROS  (20A4) 

C 

C  NOAT, IPNCHB  (214) 

C  NOATO  NUMBER  OP  ATOMS  (<«20) . 

C  IPNCHBO  PUNCH  B  MATRIX  IF  NON-ZERO  (SEE  MOTE  BELOW) . 

C 

C  X,Y,Z,ID  (3G12.6.11A4) 

C  X,y,Z0  CARTESIAN  COORDINATES  OF  AN  ATC»1. 

C  IDO  FREE  FORMAT  LABEL  (COLS  37-80) . 

C  REPEAT  NOAT  TIMES. 

C 

C  ICODE,I,J,K,L,IX,JX,FACTOR,ID  (7I4,G12.6, 10A4) 

C  ICODEO  INTERNAL  COORDINATE  TYPE  (SEE  BELOW).  IP  ICODE<0,  THE 

C  NEW  B  MATRIX  ELEMENTS  ARE  ADDED  TO  THE  PREVIOUS  ONES. 

C  I,J,K,L0  ATOM  NUMBERS  INVOLVED. 

C  IX.JXO  OPTIONAL  WEIGHTING  OF  INTERNAL  COORD  BY  THE  IX-JX  BOND 

C  LENGTH  (NOT  USED  IF  IX  AND/OR  JX  IS  0 ) . 

C  FACTORO  NEW  ROW  OF  B  IS  MULTIPLIED  BY  FACTOR  (BEFORE  BEING 

C  ADDED  TO  PREVIOUS  ROW,  IF  ICODE<0) .  FACTOR  DEFAULTS  TO 

C  1.0.  USE  TO  COMBINE  INTERNAL  COORDINATES,  IF  DESIRED. 

C  IDO  FREE  FORMAT  LABEL  (COLS  41-80) . 

C  REPEAT  AS  OFTEN  AS  REQUIRED,  TERMINATING  WITH  A  BLANK  CARD. 

C  THE  TOTAL  NUMBER  OF  INTERNAL  COORDS  MUST  BE  <=  3*NOAT. 

C 

C  ENTIRE  DECK  MAY  BE  REPEATED 
C 

C  ICODE  MODE 

C 

C  1  BOND  STRETCH 

C  I  AND  J  ARE  ATOMS  INVOLVED.  K,L,IX,JX  MUST  BE  0. 

C  2  VALENCE  ANGLE  BEND 

C  I  AND  K  ARE  TERMINAL  ATOMS,  J  IS  CENTRAL  ATOM.  L  MUST  BE  0. 

C  I,J,K  MUST  NOT  BE  COLINEAR. 

C  3  OUT  OF  PLANE  WAG 

C  I  IS  WAGGED  ATOM,  J  IS  APEX  ATC»1,  K  AND  L  ARE  ANCHOR  ATOMS. 

C  4  TORSION 

C  J  AND  K  DEFINE  THE  BOND  UNDER  TORSION.  I  AND  L  ARE  THE  NO  OF 

C  ATC»1S  (<=5)  ATTACHED  TO  J  AND  K  RESPECTIVELY.  THE  FOLLOWING  2 

C  CARDS  GIVE  THE  ATOM  NOS  FOR  THE  I-TYPE  AND  L-TYPE  ATOMS  (EACH 

C  CARD  IS  514)  .  NONE  OF  THE  I'S  OR  L'S  SHOULD  BE  THE  SAME,  OR 

C  EQUAL  TO  J  OR  K.  THE  TORSION  IS  PROPERLY  NORMALIZED  (SEE  R  L 

C  HILDERBRANDT,  J  MOLEC  SPEC,  44,  599  (1972)  )  . 

C  5  LINEAR  BEND  (DEFINES  2  INTERNAL  COORDINATES) 

C  6  LINEAR  BEND  (DEFINES  1  INTERNAL  COORDINATE) 

C  I  AND  K  ARE  END  ATOMS,  J  IS  CENTRAL  ATOM.  THE  FOLLOWING  CARD 

C  GIVES  A  POINT  (IN  3G12.6  FORMAT)  PERPENDICULAR  TO  I-J-K  AT  J 
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C  WHICH  ORIENTS  THE  BENDING  COORDINATE.  L  MUST  BE  0. 

C  FOR  ICOOE«5  A  PERPENDICULAR  INTERNAL  COORD  IS  ALSO  DEFINED. 

C 

C  B  IS  (NOB.NA)  WHERE  NA«3*NOAT 

C  B  MUST  BE  DEFINED  AS  A  SQUARE  MATRIX  FOR  PROGRAMS  2  (F  TRY/A1X»« 

C  DISP)  AND  3  (FORCE  CONSTANT  FITTER)  OF  THE  VIBRATIONAL  PACKAGE. 

C  XO  X,  Y,  Z  COORDINATES  OF  THE  ATCBIS  (SIZED  (3,NOAT)  ) 

C 

C  REQUIRED  SUBROUTINESO  BOST.  BEND,  OPLA,  TORS,  LIBE 
C 

C  SUBROUTINES  BOST,  BEND,  OPLA  AND  LIBE  WERE  MODIFIED  FROM  J  H  SCHACHT- 
C  SCHNEIDER'S  'GMAT'  PROGRAM  (SHELL  DEVELOPMENT  CO)  WITH  PERMISSION. 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

C  TO  REDIMENSION,  CHANGE  FOLLOWING  CARO  AND  ALL  OTHER  BLANK  COMMON 
COIMON  IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,70),B(200,200) 

INTEGER  TITLE(40),IDC(11),IDQ(10) 

C  READ  TITLE  CARDS 

10  READ (5. 1000, END*210 ) TITLE, NOAT,IPNCHB 
WRITE(6. 1010)TITt  NOAT 
NAaNOAT*3 
DO  20  1=1, NOAT 

READ(5,1020) (X(J, I) , J=l, 3) , IDC 
20  WRITE(6,1030)I, (X( J, I) , J=l, 3 ) , IDC 
NOB=0 
IER=0 

C  ISCAN  IS  0  NORMALLY,  >0  FOR  ERROR  SCAN  AFTER  AN  INPUT  ERROR  IS  FOUND 
ISCAN=0 
WRITE(6, 1040) 

C  READ  INTERNAL  COORD  DEFINITIONS 

3  0  READ (5,1050) ICODE , N1 , N2 , N3 , N4 , N5 , N6 , FACTOR, IDQ 
IF  (ICODE. EQ.O)  GO  TO  150 
C  IF  THIS  IS  A  NEW  COORDINATE,  INCREMENT  NOB 
IF ( ICODE , GT . 0 ) NOB=NOB+ 1 
C  DECREMENT  BY  1  IF  ICODE  =  -5 
IF ( ICODE . EQ . -5 ) NOB=NOB-l 
IF ( FACTOR . EQ . 0 . DO ) FACTOR=l . DO 

WRITE ( 6 , 1 0 6 0 ) NOB , ICODE , N1 , N2 , N3 , N4 , N5 , N6 , FACTOR , IDQ 
C  IF  THIS  ROW  IS  TO  BE  ADDED  TO  PREVIOUS  ROW,  STORE  NEW  ROW 
C  TEMPORARILY  IN  ROW  NOB+1  OF  B 
IF ( ICODE . LT . 0 ) NOB=NOB+ 1 

C  INCREMENT  BY  2  IF  ICODE=-5  SINCE  ICODE=5  DEFINED  2  ROWS  OF  B 
IF ( ICODE . EQ . -5 ) NOB=NOB+l 
IF(NOB.GT.NA)GO  TO  180 
C  ZERO  ROW  OF  B 

DO  40  J=1,NA 
40  B(NOB, J)=0.D0 
IC=IABS( ICODE) 

GO  TO  (1,2,3,4,5,6) ,IC 
WRITE(6, 1070) ICODE 
GO  TO  200 

1  CALL  BOST 
GO  TO  60 

2  CALL  BEND 
GO  TO  60 

3  CALL  OPLA 
GO  TO  60 

4  CALL  TORS 
GO  TO  60 

C  ZERO  EXTRA  ROW  OF  B  IF  ICODE=+-5 


5  IrnHOB*! 

IF(I.iST.NA)GO  TO  180 
DO  50 

50  B(I,J)«0.D0 

6  CALL  LIBE 

60  1F(IER.NE.O)GO  TO  190 
C  MULTIPLY  NEW  ROW(S)  BY  FACTOR  (IF  NOT  1.0) 

IF(FACTOR.EQ.1.DO)GO  TO  105 
IF(IC.BQ.5)GO  TO  80 

70  ISW*0 
IsNOB 
GO  TO  90 

80  ISWsl 
I«NOB-l 

90  DO  100  Jsl.NA 
100  B(I,J)=B(I,J) ‘FACTOR 
IF(ISW.EQ.l)GO  TO  70 

C  DO  WE  ADD  CURRENT  ROW(S)  TO  PREVIOUS  ROW(S)  ? 

105  IF(ICODE.GT.0)GO  TO  30 
IF(ICODE.EQ.-5)GO  TO  110 
ISWsO 
IsNOB-1 
GO  TO  130 
110  ISW=1 
120  IsNOB-2 
130  DO  140  J=1,NA 
140  B(I, J)=B(I,J)  +  B(NOB,J) 

NOB=NOB-l 

IF{ISW.EQ.0)GO  TO  30 

ISW«0 

GO  TO  120 

150  IF(lSCAN.NE.0)GO  TO  10 
WRITE (6, 1080) NOB 
K»-ll 
160  K=K-^12 

L=MIN0(K+11,NA) 

WRITE(6, 1090) (J, J=K,L) 

DO  170  1=1, NOB 

OPEN  { 22 ,  FILE= '  BMAT .  IN ' ) 

WRITE(22,1101)  (B(l, J) . J=K,L) 

170  WRITE(6, 1100)1, (B{I,J),J=K,L) 

C. .DZ  DO  171  1=1,15 

C..DZ  WRITE(22,1101)  (B(l, J) , J=l, 12) 

C..DZ  171  CONTINUE 
C..DZ  DO  172  1=1,15 
C..DZ  WRITE(22,1101)  (B(I, J) , J=13,21) 

C..DZ  172  CONTINUE 

1F(L.LT.NA)G0  TO  160 

C  EACH  ELEMENT  OF  B  IS  PUNCHED  IN  A8  FORMAT  -  THE  INTERNAL  64  BIT  (8 
C  BYTE)  FLOATING  POINT  NUMBER  IS  INTERPRETED  AS  8  EBCDIC  CHARACTERS  (1 
C  CHARACTER  IS  STORED  IN  1  BYTE  (=  8  BITS)  IN  IBM  360/370  COMPUTERS)  . 

C  EACH  DOUBLE  PRECISION  (REAL‘8)  VALUE  THEN  OCCUPIES  8  CARD  COLUMNS  - 
C  THIS  FORMAT  MINIMIZES  THE  SIZE  OF  THE  B  MATRIX  CARD  DECK,  BUT  IS 
C  THEN  COMPLETELY  INCOMPREHENSIBLE.  DO  NOT  INTERPRET  THESE  CARDS. 

IF(IPNCHB.NE.0)WRITE(7,1160)TITLE,NOB.NA, ( (B(I, J)  .  I=l,NOB) , J=1,NA) 
GO  TO  10 

180  WRITE(6,1140} 

STOP 

190  1F(IER.EQ.1)WRITE(6,1130) 

200  IF(ISCAN.EQ.0)WRITE(6,1170) 
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C  ERROR  SCAM  FOR  THIS  DATA  DECK,  AND  DON'T  PRINT/PUMCH  B  MATRIX 
ISCANxISCAN-t-l 
lERsO 
GO  TO  30 

210  WRITE(6,1120) 

STOP 

1000  FORMAT(20A4/20A4/214) 

1010  FORMAT('1',20A4,24X, 'BMAT  (VERSIONO  JVL  28,  1977) ' /IX, 20A4/ 

$  'ONUMBER  OP  ATOMS  , 14/ ' OATCW' , 8X, 'X' , IIX, 'V' , IIX, 'Z' , IIX, 'ID' ) 

1020  FORMAT (3G12. 6, 11A4) 

1030  FORMAT( ' 0 ' , 13 , 2X, 3F12 . 6, 5X, 11A4) 

1040  FORKAT(/'0  INTERNAL  COORDINATE  DEFINITIONSO ' / ' ONOB  CODE  I 
$  'J  K  L  IX  JX  FACTOR') 

1050  FORMAT(7I4,G12.6,10A4) 

1060  FORMATdX,  I3,2I5,5I4,F11.6,1X,10A4) 

1070  FORMAT{ ' OILLEGAL  CODE' , IS, '  CHOSEN' ) 

1080  FORMAT ('ONUMBER  OF  INTERNAL  COORDINATES  s', 14/ 'IB  MATRIX  ', 

$  '(NOB  BY  3*NOAT)0') 

1090  FORMAT(//3X, 12110) 

1100  FORMAT('0',I4,2X,12F10.6) 

CllOl  FORMAT(I4,2X, 12E15.6) 

1101  FORMAT(12E15.6) 

1120  FORMATCl***  NORMAL  TERMINATION'//) 

1130  FORMAT('  ILLEGAL  SPECIFICATION  OF  I,  J,  K,  L,  IX  OR  JX') 

1140  FORMATCO***  PROGRAM  TERMINATED  -  TOO  MANY  INTERNAL  COORDS'//) 
C1160  FORMAT(20A4/20A4/2I4/(10A8)) 

1160  FORMAT(20A4/20A4/2I4/ (3D23 .16) ) 

1170  FORMAT('0***  PROGRAM  WILL  SCAN  FOR  FURTHER  ERRORS  IN  DATA  DECK'/) 
END 

SUBROUTINE  BOST 

C  ‘THIS  SUBROUTINE  COMPUTES  THE  B  MATRIX  ELEMENTS  FOR  A  BOND  S'^RETCH 
C  AS  DEFINED  BY  WILSON. 

IMPLICIT  REAL*8  (A-H,0-2) 

COMMON  IC, I, J,K,L, IX, JX, NOAT,NOB, IER,X(3, 70) , B (200, 200) 

COMMON/SCHACH/RIJ(3) ,RJK(3) , RJL(3) , EIJ(3) , EJK(3 ) , EKL(3 ) 

IP(I.LE.O.OR.I.GT.NOAT)GO  TO  30 

IF(J.LE.O.OR.J.GT.NOAT)GO  TO  30 

IF(K.NE.0)GO  TO  30 

IF(L.NE.0)GO  TO  30 

IF(IX.NE.0)GO  TO  30 

IF(JX.NE.0)GO  TO  30 

DIJSQsO.DO 

DO  10  M=l,3 

T=X(M,J)-X(M, I) 

RIJ(M)=T 

10  DIJSQeDIJSQ+T*T 
DlJsDSQRT(DIJSQ) 

11*3* (I-l) 

JJ=3*(J-1) 

DO  20  M=1.3 
T*RIJ(M) 

IF(DABS(T) .LT.l.D-8)GO  TO  20 
TsT/DlJ 

B(NOB,II+M)=-T 
B(NOB, JJ+M)=T 
20  CONTINUE 
RETURN 
30  lERsl 
RETURN 
END 
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SUBROUTINE  BEND 

C  THIS  SUBROUTINE  OXIPUTES  THE  B  MATRIX  ELQ<ENTS  OF  A  VALENCE 

C  ANGLE  BENDING  COORDINATE  AS  DEFINED  BY  WILSON. 

C  I  AND  K  ARE  THE  NUMBERS  OF  THE  END  ATOMS. 

C  J  IS  THE  NUMBER  OF  THE  CENTRAL  ATC»<. 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  IC, I , J, K, L, IX. JX. NOAT. NOB. lER. X{3,70),B(200,200) 

CCMMON/SCHACH/RJI ( 3 ) . RJK ( 3 ) . RJL ( 3 ) . EJI ( 3 ) . EJK ( 3 ) . EKL ( 3 ) 

IF(I.LE.O.OR.I.GT.NOAT}GO  TO  50 

IP(J.LE.O.OR.J.GT.NOAT)GO  TO  SO 

IP(K.LE.O.OR.K.GT.NOAT)aO  TO  50 

IF(L.NE.O)GO  TO  SO 

IF(IX.LT.O.OR.IX.GT.NOAT)GO  TO  50 

IF(JX.LT.O.OR.JX.GT.NOAT)GO  TO  50 

IF(IX.NE.O.AND.JX.NE.O)GO  TO  10 

IXsl 

JXsl 

10  DJISQsO.DO 
DJKSQsO.DO 
DXSQ=0.D0 
DO  20  M=l,3 
TP=X(M. J) 

TaX(M.I)-TP 

RJI (M) =T 

DJISQ=DJISQ+T*T 

TaX(M,K)-TP 

RJK (M) =T 

DJKSQ=DJKSQ+T*T 

T*X(M,JX)-X(M,IX) 

20  DXSQ*DXSQ+T*T 
DCriaDSQRT(DJlSQ) 

DJK»DSQRT(DJKSQ) 

DX=DSQRT(DXSQ) 

IF(DX.EQ.0.D0)DX=1.D0 
DOTJsO.DO 
DO  30  M=l,3 
TsRJI (M) /DJI 
EJI(M) =T 
TP=RJK(M)/DJK 
EJK(M)=TP 
30  DOTJ=DOTJ+T*TP 

IF(DABS(DOTJ) .GT.0.99995D0)GO  TO  60 
SINJsDSQRT ( 1 .DO-DOTJ*DOTJ) 

II=3*(I-1) 

JJ*3*(J-.l) 

KK=3*(K-1) 

DO  40  Msl.3 

SMI*DX*(DOTJ*EJI (M) -EJK(M) ) / (DJ1*SINJ) 

IF(DABS(SMI)  .GE.l.D-8)B(NOB.II-i-M)<SMI 
SMK*DX* {DOTJ*EJK(M) -EJI (M) ) / (DJK*SINJ) 

IF{DABS(SMK) .GE. 1 .D-8) B (NOB.KK+H) sSMK 
SUMsSMI-fSMK 

40  IF(DABS(SUM} .GE. 1 .D-8) B(NOB. JJ+H) s-SUM 
RETURN 
50  IER=1 
RETURN 
60  IER=-1 

WRITE (6, 1000) 

RETURN 

1000  FORMATC  I-J-K  IS  COLINEAR  -  USE  LINEAR  BEND') 


END 

SUBROUTINE  OPLA 

C  THIS  SUBROUTINE  COMPUTES  THE  B  MATRIX  ELQ4ENTS  FOR  AN  OUT  OF 
C  PLANE  WAGGING  COORDINATE  AS  DEFINED  BY  DECIUS.  MCINTOSH,  MICHAELIAN 
C  AND  PETERSON.  SUBROUTINE  CODED  BY  M  PETERSON,  UNIV  OF  TORONTO. 

C  I  IS  THE  END  ATOM  (ATC»!  WAGGED  WITH  RESPECT  TO  J-K-L  PLANE)  . 

C  J  IS  THE  APEX  ATOM  (ATOMS  I,  K  AND  L  ARE  ATTACHED  TO  J) . 

C  K  AND  L  ARE  THE  ANCHOR  ATOMS  (DEFINE  THE  J-K-L  PLANE) . 

IMPLICIT  REALMS  (A-H.O-Z) 

COMMON  IC.  I, J,K.L, IX. JX,NOAT.MOB. IER,X(3, 70) , B(200, 200) 
COMMON/SCHACH/RJI ( 3 ) , RJK ( 3 ) , RJL ( 3 ) . EJI ( 3 ) ,  EJK ( 3 ) , EJL ( 3 ) 

DIMENSION  Cl (3) 

IP(I.LE.0.OR.I.GT.NOAT)GO  TO  60 
IP(J.LE.O.OR.J.GT.NOAT)GO  TO  60 
IF(K.LE.O.OR.K.GT.NOAT)GO  TO  60 
IF(L.LE.O.OR.L.GT.NOAT)GO  TO  60 
IP(IX.LT.O.OR.IX.GT.NOAT)GO  TO  60 
IF(JX.LT.O.OR.JX.GT.NOAT)GO  TO  60 
IF(IX.NE.0.AND.JX.NE.0)GO  TO  10 
IXal 
JX=1 

10  DJISQ=0.D0 
DJKSQ=0.D0 
DJLSQsO.DO 
DXSQaO.DO 
DO  20  M=l,3 
TP=X(M, J) 

TsX(M,I)-TP 

RJI (M) =T 

DJISQ=DJISQ+T*T 

T*X(M,K)-TP 

RJK(M)=T 

DJKSQ=DJKSQ+T*T 

T=X{M,L)-X(M, J) 

RJL(M)=T 

DJLSQ=DJLSQ+T*T 

T*X(M,JX)-X(M,IX) 

20  DXSQ=DXSQ+T*T 
DJI=DSQRT(DJISQ) 

DJK=DSQRT(DJKSQ) 

DJL=DSQRT(DJLSQ) 

DX=DSQRT(DXSQ) 

IF(DX.EQ.0.D0)DX=1.D0 
COSIsO.DO 
COSK*O.DO 
COSLaO.DO 
DO  30  M=l,3 
TsRJI (M) /DJI 
EJI (M) =T 
TP=RJK(M) /DJK 
EJK (M) =TP 
TPP=RJL(M)/DJL 
EJL (M) =TPP 
COSI=COSI+TP*TPP 
COSK=COSK+T*TPP 
30  COSLsCOSL+T*TP 

IF(DABS(COSI) .GT.0.99995D0)GO  TO  70 
S INS IN= 1 . DO -COS I *COS I 
SINIaDSQRT(SINSIN) 

C1(1)5:EJK(2)*EJL(3)-EJK{3)*EJL(2) 
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Cl (2 ) «EJK ( 3 ) *EJL ( 1 ) -EJK ( 1 ) *EJL(3 ) 

C1(3)«BJK(1)*EJL(2)-EJK(2)*EJL(1) 

DOT.EJI { 1 ) *C1 ( 1 ) ♦EJl ( 2 ) -Cl ( 2 ) ♦EJI ( 3 ) *C1 { 3 ) 

SINT«DOT/SINI 

IF(DABS(SINT) .GT. 0 . OOOOSDO)WRZTE(6, 1020) 

IF(DXBS(SINT) .GT.0.99995D0)GO  TO  80 
COST*DSQRT ( 1 . D0-SINT*SINT) 

TANT»SINT/COST 

II«3*(I-1) 

JJ»3*(J-1) 

KK*3*(K-1) 

LLs3*(L-l} 

COSSIN»COST*S INI 

DO  50  M>1,3 

T«C1 (M) /COSSIN 

SMI» (T-TANT*EJI (M) ) /DJI 

IF(DABS(SMI)  .GE.l.D-8)B(NOB,II-»^M)«DX*^I 

SMKsT* (COSI*COSK-COSL) / (SINSIN*DJK) 

IF  (DABS  (SMK)  .GE .  1  .D-8 )  B  (NOB,  KK-i^M)  «DX*SMK 
SMLsT* (COSI*COSL-COSK) / (SINSIN*DJL) 

IF(DABS(SML)  .GE.l.D-8)B(NOB,LL-».M)sDX*SML 
SUMsSMI-fSMK-fSML 

50  IF(DABS(SUM)  .GE.  1  .D-8) B (NOB,  JJ-fM) x-DX*SUM 
RETURN 
60  lERsl 
RETURN 
70  IER=-1 

WRITE (6, 1000) 

RETURN 
80  IER*-1 

WRITE(6,1010) 

RETURN 

1000  FORMAT ('  K-J-L  IS  COLINEAR  (NO  FLANE  DEFINED  FOR  WAG  OF  I)') 

1010  FORMAT ('  I  IS  PERPENDICULAR  TO  J-K-L  PLANE  -  USE  VALENCE  ANGLE  ', 

$  'BENDS') 

1020  FOBMAT('+',86X, '**•  WARNINGO  WAG  OF  A  NON-PLANAR  SYSTEM  ***') 

END 

SUBROUTINE  TORS 

C  THIS  SUBROUTINE  COMPUTES  THE  B  MATRIX  ELEMENTS  FOR  TORSION  AS  DEFINED 
C  BY  R  L  HILDERBRANDT  IN  J  MOLEC  SPEC,  44,  599  (1972) . 

C  SUBROUTINE  CODED  BY  M  PETERSON,  DEPT  OF  CHQtlSTRY,  UNIV  OF  TORONTO. 

C  J  AND  K  DEFINE  THE  BOND  UNDER  TORSION. 

C  NI  AND  NL  ARE  THE  NUMBER  OF  ATOMS  ATTACHED  TO  J  AND  K  RESPECTIVELY 
C  (NI,  NL  <=  5) .  2  DATA  CARDS  ARE  READO  (1)  CONTAINS  NI  ATOM  NUMBERS 
C  FOR  THE  I-TYPE  ATOMS,  AND  (2)  CONTAINS  NL  ATOM  NUMBERS  FOR  THE  L-TYPE 
C  ATOMS  (BOTH  CARDS  ARE  IN  514  FORMAT) . 

C 

C  lATOM,  LATOMO  ATOM  NUMBERS  FOR  THE  I-  AND  L-TYPE  ATOMS  (SIZED  5) 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  IC,NI, J,K,NL,IX, JX,NOAT,NOB,IER,X(3,70) ,B(200,200) 
CC»1MON/SCHACH/RI J ( 3 ) , RJK ( 3 ) , RLK ( 3 ) . El J (3 ) , EJK ( 3 ) , ELK ( 3 ) 

DIMENSION  CR(3),IATOM(5) ,LATOM(5),SJ(3),SK(3) 

READ(5, 1000) (lATOM(I) ,1*1, NI) 

WRITE{6,1010) (lATOM(I) ,I*1,NI) 

READ(5,1000) (LATOM(L),L=l,NL) 

WRITE(6,1020) (LATOM(L) ,L=1,NL) 

IF(NI.LE.O.OR.NI.GT.5)GO  TO  110 
IF(J.LE.O.OR.J.GT.NOAT)GO  TO  110 
IF(K.LE.O.OR.K.GT.NOAT)GO  TO  110 
IF(NL.LE.O.OR.NL.GT.5)GO  TO  110 
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XF(XX.LT.O.OR.IX.GT.HOAT)GO  TO  110 
XF(JX.LT.O.OR.JX.GT.NOAT)GO  TO  110 
XF(XX.NE.O.AND.JX.NE.O)GO  TO  10 
XX«1 
JX«1 

10  DJKSQsO.DO 
DXSQ«0.D0 
DO  20  Msl,3 
SJ<M)*0.D0 
SK(M) «0.D0 
T«X(M,K)-X(M,J) 

RJK(M} «T 

DJKSQ»DJKSQ+T*T 

T>X(M.JX)-X(M,XX) 

20  DXSQ«DXSQ>T*T 

OJKsl .O0/DSQRT(OJKSQ) 

DX«DSQRT(DXSQ) 

XF (DX . EQ . 0 . DO } DXsl .DO 
DO  30  Msl.3 

30  EJK{M)*RJK(M)*DJK 
JJa3*(J-l) 

KK=3*(K-1) 

C  LOOP  OVER  THE  X-TYPE  ATOMS 
DO  60  Nsl.NI 
XsXAT(»l(N) 

XF(I.LE.O.OR.X.GT.NOAT)GO  TO  110 
DXJSQsO.DO 
DO  40  Msl,3 
T»X(M, J)-X(M, X) 

RIJ(M)sT 

40  DIJSQ»DIJSQ+T*T 

DXJsl . DO /DSQRT (DIJSQ) 

COSJsO.DO 
DO  50  Msl,3 
T=RXJ(M)*DIJ 
EXJ(M)=T 

50  COSJ=COSJ-T*EJK(M) 

XF(DABS(COSJ) .GT.0.99995D0)GO  TO  120 
SIN2J=(1.DO-COSJ*COSJ) *DFLOAT(NI) 
XI=3*(X-1) 

CR(1)=EXJ{2)*EJK{3)-EXJ(3)*EJK(2) 

CR(2)=EXJ{3)*EJK(1)-EXJ{1)*EJK{3) 

CR{3)*EXJ(1)*EJK(2)-EXJ{2)*EJK(1) 

DO  60  Msl,3 
TaCR(M) /SIN2J 
SMX*T*DXJ 

XF  (DABS  ( SMX )  .GE .  1 . D-8 )  B  (NOB,  XX-t-M)  s-DX*SMI 

aa*T*COSJ*DJK 

SK(M)sSK(M)4-SMK 

SMJsSMX-SMK 

60  SJ  (M)  »SJ  (M) -fSMJ 
C  LOOP  OVER  THE  L-TYPE  ATOMS 
DO  90  N=1,NL 
L>LATOM(N) 

XF(L.LE.O.OR.L.GT.NOAT)GO  TO  110 
DLKSQsO.DO 
DO  70  M=l,3 
T*X(M,K)-X(M,L) 

RLK(M)sT 

70  DLKSQaDLKSQ+T*T 
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DLK-1 . OO/OSQRTtOLKSQ) 

COSK>O.DO 
DO  80  M>1,3 
T«RLK(M) *DLK 
ELK(M)«T 

80  COSK-COSX^EJK(M}*T 

IP(DRBS(COSK) .GT.0.99995D0)GO  TO  120 
SIN2K> ( 1 . DO-COSK*COSK) *OFU>AT (NL) 

LL-3*(L-1) 

CR ( 1 ) -ELK ( 3 ) * EJK ( 2 ) -ELK ( 2 ) * EJK ( 3 ) 
CR(2)«ELK(1)*EJK(3)-ELK(3)*EJK(1) 

CR ( 3 ) «ELK (2 ) *EJK ( 1 ) -ELK ( 1 ) *EJK ( 2 ) 

DO  90  M«l,3 
T«CR(M) /SIN2K 
SML«T*DLK 

ZF(DABS(SML) .GE.l.D-8)B(NOB,LL^M)*-DX*SML 
SMJ«T*COSK*DJK 
SJ(M}sSJ(M)-»^SHJ 
SMKsSML-SMJ 
90  SK(M)sSK(M)4-SMK 
DO  100  M»l,3 
SMJsSJ  (M) 

IF  ( DABS  ( SMJ )  . GE .  1 . 0-8 )  B  (NOB .  JJ-i-M)  «SMJ*DX 
SMK«SK(M) 

100  IF(DABS(SKK)  .GE.  1  .D-8)  B(NOB, KK-*-H)  «SMK*DX 
RETURN 
110  IER=1 
RETURN 
120  IER*-1 

WRITE(6,1030) 

RETURN 

1000  F0RMAT(SI4) 

1010  FOFMAT{'+',86X, 'I0',5I4) 

1020  FORMAT('  +  M09X, ',L0',5I4) 

1030  FORMAT ('  I-J-K  OR  J-K-L  IS  COLINEAR  (NO  TORSION  POSSIBLE)') 
END 

SUBROUTINE  LIBE 

C  THIS  SUBROUTINE  COMPUTES  THE  B  MATRIX  ELB1ENTS  FOR  A  LINEAR  BEND 
C  OR  FOR  A  PAIR  OF  PERPENDICULAR  LINEAR  BENDS. 

C  I  AMD  K  ARE  THE  END  ATOMS. 

C  J  IS  THE  CENTRAL  ATOM. 

C 

C  A  GIVES  THE  CARTESIAN  COORDINATES  OF  A  POINT  IN  SPACE,  SUCH 
C  THAT  THE  VECTOR  FROM  AT(»i  J  TO  POINT  A  IS  PERPENDICULAR  TO 
C  THE  LINE  I-J-K  AMD  SERVES  TO  ORIENT  THE  COORDINATES  IN  SPACE. 
IMPLICIT  REAL*8  (A-H,0-Z) 

CCMMON  IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,70)  ,B(200,200) 
CC»IMON/SCHACH/RJI  ( 3 ) ,  RJK  ( 3 ) ,  EJK  ( 3 ) ,  UP  ( 3 ) ,  UN  ( 3 ) ,  UNIT  ( 3 ) 
DIMENSION  A(3) 

READ(5,1000)A 
WRITE (6. 1010)  A 

IF(I.LE.O.OR.I.GT.NOAT)GO  TO  SO 
1F(J.LE.O.OR.J.GT.NOAT)GO  TO  60 
IF(K.LE.O.OR.K.GT.NOAT)GO  TO  60 
1F(L.NE.O)GO  TO  60 
IF(IX.LT.O.OR.1X.GT.NOAT)GO  TO  60 
1F(JX.LT.O.OR.JX,GT.NOAT)GO  TO  60 
1F(IX.NE.O.AND.JX.NE.O)GO  TO  10 
IXsl 
JXsl 
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10  DJISQ«0.D0 
DJKSQ«O.DO 
DXSQ«0.D0 
DJASQ«0.D0 
DO  20  M«1.3 
TP«X(M,J) 

T«X{M,1)-TP 

RJI{M)»T 

DJISQ«DJISQ+T*T 

T«X(M,K)-TP 

RJK(M)«T 

DJKSQ»DJKSQ+T*T 

T«X(M,JX)-X(M,IX) 

DXSQ«DXSQ-^T*T 

T«A(M)-TP 

UN(M)«T 

20  DJASQsDJASQ4-T*T 
DJI>DSQRT(DJISQ) 

DJKsDSQRT ( DJKSQ ) 

DXsDSQRT(DXSQ) 

DJAsDSQRT(DJASQ) 

IF(DX.EQ.0.D0)DXsl.D0 

DOTJsO.DO 

DOTP«O.DO 

DO  30  Msl,3 

TssRJI  (M)  /DJI 

TP*RJK(M)/DJK 

EJK(M)=TP 

DOTJ=DOTJ+T*TP 

TP*UN(M)/DJA 

UNIT(M)sTP 

30  DOTP=DOTP+T*TP 

TESTsDABS ( DOTJ ) -1 . DO 
IF(DABS(TEST) .GT.0.00005D0)GO  TO  70 
1F(DABS(D0TP) .GT.0.00005D0)GO  TO  80 
II=3*{I-1) 

JJ*3*(J-1) 

KK»3*(K-1) 

DO  40  M=l,3 
T=UNIT(M) 

IF(DABS{T) .LT.l.D-8)GO  TO  40 

T=-DX*T 

SMI=T/DJI 

6(NOB,II-»-M)sSMI 

SMKsT/DJK 

B(NOB.KK-i-M)«SMK 

B  (NOB ,  JJ-i-M)  s-^I-SMK 

40  CONTINUE 

IF(IC.EQ.6)RETURN 

NOBsNOB-fl 

UP { 1 ) =EJK (2 ) *UNIT ( 3 ) -EJK { 3 ) *UNIT{2) 
UP ( 2 ) *EJK ( 3 ) *UNIT ( 1 ) -EJK ( 1 ) *UNIT ( 3 ) 
UP ( 3 ) «EJK ( 1 ) *UNIT ( 2 ) -EJK ( 2 ) *UNIT ( 1 ) 
DO  50  Msl,3 
T*DP(M) 

IF(DABS(T) .LT.l.D-8)GO  TO  50 

T*-DX*T 

SMI*T/DJI 

B  (tK)B,  II-t-M)  cSMI 

SMKsT/DJK 
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B(NOB.KK-fM)«SMK 
B  (NOB.  JJ4-M)  «-SMI-SMK 
SO  CONTINUE 
RETURN 
60  ZER«1 
RETURN 
70  IER«-1 

WRITE (6, 1020} 

RETURN 
80  IER«>1 

WRITE(6.1030) 

RETURN 

1000  FORHRT(3G12.6) 

1010  FORMAT('  +  ',86X, 'A  =  ( ' .2 (F11.7. ' , , F11.7, ' ) ' ) 

1020  FORKATC'  I-J-K  NOT  COLINEAR  >  USE  VALENCE  ANGLE  BEND') 
1030  FORMAT ('  ATOM  A  NOT  PERPENDICULAR  TO  I-J-K  AT  J') 

END 
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IL2.2  Program  matmultf 

The  next  program  is  matmultf.  This  program  carries  out  the  transformatitHi 

F=B-*KB-‘  (7) 

and  also  die  detaminadon  of  K  from  F.  The  parameter  NAT  in  the  program  rqnesents  the 
number  of  attmis  in  the  molecule  considered  and  must  be  changed  for  each  molecule  considered. 
The  file  containing  the  matrix  K.  KMATJN,  can  be  obtained  from  a  Gaussian  calculation  or  a 
C^PAC  calculation. 

On  the  following  pages  a  listing  of  the  FORTRAN  program  matmultf  is  given. 
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matault . f 


PROGRAM  MAIN 

PARAMETER (NATsl6 , MMs3 •HAT-6, Ns3*NAT, MMM>2 •MM, NROWsMM, 

♦  NMATR=NROW^ (NROW* 1 ) /2 ) 

IMPLICIT  REALMS  (A-H,0-Z) 

REALMS  K(N.N) ,M(N,N) .FID(NMATR) 

REAL  MC,MO.MH 

DIMENSIOH  B(MM,N)  ,BM(MM.N)  ,BP(N,M<)  ,F(MM,MH)  ,TBST(MM,M<) , 

1  BMBP  (MM,  MM) .  BMBPI  (MM,  MM) ,  BPIM  (Ml,  N) ,  PROD  (MM,  N) .  BPLMI  (N,  Ml) . 

2  TESTA  (MM,  MM) ,  AA  (MM,  NMM) ,  BB  (Ml,  MM)  ,  TEST2  (MM,  MM) 

CCMMON  NSYS. INDEX, DET 

C. .  . 

C .  . .  GET  MAlltlCES  B  AND  K 
C. .  . 

CALL  BKMATR(MM,N,B,K) 

C. .  . 

C . . ,  ADJOINT  OF  B 
C.  .  . 

DO  10  1*1, N 
DO  10  Jsl.MM 

BP(I,J)=B(J.I) 

10  CONTINUE 
C.  ,  . 

C . . .  DETERMINE  PRODUCT  OP  B  M  BP  MATRICES 
C.  .  . 

MC*12.01 
MO*16.00 
MHsl.OOS 
DO  501  1*1, N 
DO  SOI  J*1,N 
M(I, J)*0. 

501  CONTINUE 

DO  502  1=1, N 
M(I,I)=1. 

502  CONTINUE 
C.  .  . 

DO  11  1*1, MM 
DO  11  J=1,N 

SUM=0. 

DO  11  L=1,N 

SUM*SUM+B ( I , L) *M (L, J) 

BM(I, J)*SUM 

11  CONTINUE 
C.  . . 

C. . .  DETERMINE  PRODUCT  OF  B  M  BP  MATRICES 

DO  511  1*1, MM 
DO  511  J*1,MM 
SUM*0. 

DO  511  L*1,N 

SUM*SUM-^BM  ( I ,  L)  •BP  (L,  J) 

BMBPd,  J)sSUM 

511  CONTINUE 

OPEN ( 3 , FILE* ' TESTADZ . OUT ' ) 

OPEN ( 23 , FILE* ' BMBP .mat ' ) 

WRITE(3,^)  'BMBP' 

DO  540  1*1, MM 

WRITE(3,115)  (BMBPd,  J),J*1,MM) 

WRITE(23,115)  (BMBPd,  J)  ,J»1, MM) 
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540  CC»<T1NUE 

C.  . . 

C. . .  DETERMINE  INVERSE  OF  B  M  BP 

DO  12  I«1.MM 
DO  12  J«1,MM 

AAd.  J)«BMBP(I,J) 

12  CONTINUE 
NSYS«0 
INDEX*! 

DO  221  I«l,)&f 
DO  221  J«MM-fl,>MM 
AA(I,J)«0. 

IF((J-M<) .EQ.I)  AA(1.J)*1. 

221  CONTINUE 

DO  222  1*1, MM 
DO  222  J*1,MM 
BB(I,J)*0. 

IF(I.EQ.J)  BB(I,J)*1. 

222  CONTINUE 

CALL  MATCALC(AA.BB,MM,MMM) 

WRITE(6,*)  'DETA  a',DET 

C. .  . 

C . . .  SET  BMBPI  MATRIX 

C.  .  . 

DO  191  1*1, MM 
DO  191  J*1,MM 

BMBPI  ( I ,  J)  *AA( I,  J-fMM) 

191  CONTINUE 
WRITEO,*)  'BMBPI' 

DO  192  1*1, MM 

WRITEO, 116)  (BMBPKI,  J) ,  J*1,MM) 

192  CONTINUE 
C.  .  . 

C . . .  DETERMINE  TESTA  MATRIX 

C.  .  . 

DO  302  1*1, MM 
DO  302  J*1,MM 
SUM*0. 

DO  302  L*1,MM 

SUM*SUM-i-BMBPI  ( I ,  L)  *BMBP  (L,  J) 
TESTAd,  J)=SUM 

302  CONTINUE 

WRITE(3,*)  'TESTA' 

DO  340  1*1, MM 

WRITE(3,111)  (TESTAd,  J),J*1, MM) 

340  CONTINUE 

C  WRITE (3, 102) 

C  DO  341  1*1, MM 

C  WRITE(3,112)  (TESTAd,  J),J*13,N) 

C  341  CONTINUE 
C.  .  . 

C .  .  .  DETERMINE  BPLM  MATRIX 

C.  .  . 

DO  13  1*1, MM 
DO  13  J=1,N 

SUM*0. 

DO  13  L=1,MM 

SUM*SUM-«- BMBPI  (I ,  L)  *BM(L,  J) 
BPUfd,  J)*SUM 

13  CONTINUE 
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c. . 
c. . 


c. .  . 


202 


240 

C 

c 

c 

C  241 
C  •  •  • 

c. .  . 
c. .  . 


14 


c. .  . 


43 


44 


C.  . . 


20 
C. . 
C. , 
C. . 
C. . 


I^TBRMINE  BPLM  *  BP  MATRIX 

DO  202  Z«1,M4 
DO  202  J*1.MM 
SDM«0. 

DO  202  L«1,N 

SUM«SUM-«-BPLM ( I .  L)  *BP  ( L,  J ) 

TEST(I,J)«SUM 

CONTINUE 

OPEN ( 2 , PILE* ' TESTDZ . OUT ' ) 

WRITE(2,*)  'BPUi  *  BP* 

DO  240  1*1, S 

WRITB(2.111)  (TEST(I,J).J*1.S) 

CONTINUE 
WRITE (2. 102) 

DO  241  1*1,15 

WRITE(2,112)  (TEST(I.J),J«13,15) 

CONTINUE 

DETERMINE  TRANSPOSE  OF  BPLM  MATRIX 

DO  14  1*1, N 
DO  14  J*1,MM 

BPLMI ( I , J) «BPLM ( J, I ) 

CONTINUE 

WRITE(2,*)  *K  MATRIX 
DO  43  1*1, IS 

WRITE(2,111)  (K(I,J), Jal.l2) 

CONTINUE 
WRITE (2, 102) 

DO  44  1*1,15 

WRITE{2,112)  (K(I,J),Jsl3,15) 

CONTINUE 

DO  20  1*1, MM 
DO  20  J*1,N 

SUM*0. 

DO  20  L*1,N 

SUM*SUM-f  BPLM ( I ,  L)  *K (L,  J) 

PRODd,  J)*SUM 

CONTINUE 

DETERMINE  PRODUCT  OF  BPUf  K  BPUfI  MATRICES  TO  GIVE 
F(I,J)  MATRIX  -  UNITS  OF  HARTREES/BOHR**2 

DO  30  1*1, MM 
DO  30  J*1,MM 
SUM*0. 

DO  30  L*1,N 

SUM*SUM4PR0D ( I ,  L)  *BPLMI  (L,  J) 

F(I,J)*SUM 

TO  PUT  F(I,J)  IN  UNITS  OF  MDYNE/A 
INSERT  THE  FOLLOWING  STATEMENT 

F{I,J)al5.57*F{I, J) 


CONTINUE 
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c. .  . 

C. .  .  F(I,J)  MATRIX 
C.  .  . 

OPEN ( 1 , FILE= ' FORCE . OUT ' ) 

OPEN (31, FILE* ' FORCE2 . OUT ' ) 

WRITE (1,*)  '  FORCE  CONSTANT  MATRIX' 

WRITE (1.*)  '  (INTERNAL  COORDINATES  -  UNITS  OF  HARTREES/BOHR**2) ' 
WRITEd,*)  ' 

II>1 

DO  36  Isl.NROW 
DO  36  J«1,I 
P1D(II)*F(I,J) 

II«II-»^1 

36  CONTINUE 

CALL  OUTPAK(FlD.NROW,NMATR.l,l) 

C  DO  40  Isl.Ml 

C  WRITEd,  111)  (F(I,  J) ,  J«1,MM) 

C  40  CONTINUE 

WRITEd,  102) 

C  DO  41  1x1,15 

C  WRITEd, 112)  (F(1,J).J=13,15) 

C  41  CONTINUE 

C.  .  . 

C . . .  CONVERT  ELEMENTS  OF  FORCE  CONSTANT  MATRIX 
C . . .  TO  UNITS  OF  MDYNES/ANGSTRCai 

C.  ,  . 

DO  42  1=1, MH 

DO  42  J=1,I 

C  F(I, J)=15.56923*F(*, J) 

WRITE(31,120)  I,J,F(I,J) 

42  CONTINUE 

C  •  • . 

C . . .  FORMATS 

C.  .  . 

102  FORMAT (IX) 

111  FORMAT (12F12. 6) 

112  FORMAT(9F12.6) 

115  FORMAT (15F12.S) 

116  FORMAT (15E15. 6) 

120  FORMAT (214,020. 12) 

STOP 

END 

SUBROUTINE  BKMATR(M,N, L, K) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REALMS  K 

DIMENSION  B(M,N) ,K(N,N} 

NATxN/3 
MMx3*NAT-6 
C. . .  B  MATRIX 

OPEN  ( 21 ,  FILEx ' BMAT .IN') 

KKx-11 

160  KKxKK+12 

L=MIN0(KK+11,N) 

DO  170  Ixl.MM 

170  READ(21,114)  (B(I, J) , JxKK,L) 

IF(L.LT.N)GO  TO  160 
C  DO  30  Ixl,MM 

C  READ(21,114)  {B(I.J),Jxl,12) 

C  30  CONTINUE 
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r>  o  o 


DO  31  Ib1,MM 

READ(21.114)  (B(I.J).J*13,N) 

31  CONTINUE 

OPEN(22,FILE='BHAT.OUT' ) 

KK*-11 

161  KKsKK4l2 

L«MINO(KK-t-ll,N) 

DO  171  1*1, MM 

171  WRITE(22,114)  (B(I. J) , J*KK.L) 

IP(L.LT.N)GO  TO  161 
DO  40  1*1, MM 

WRITE(22,110)  {B(I.J),J*1.12) 

CONTINUE 
WRITE (22. 102) 

DO  41  1*1. MM 

WRITE(22.110)  (B(I.J),J»13.N) 

CONTINUE 

K  MATRIX 

OPEN (11. FILE* ' KMAT .IN') 

DO  50  1*1. NAT 

READdl.  121)  (K(J,  1*3-2)  ,K(J.  1*3-1)  ,K(J.  1*3) .  J*1,N) 
CONTINUE 
READdl. 105) 

DO  51  1*1,21 

READdl, 101)  (Kd.J),J*10,18) 

CONTINUE 
READdl,  105) 

DO  52  1*1,21 

READdl, 104)  (K(I,J).J*19,21) 

CONTINUE 

OPEN (12 , FILE* ' KMAT . OUT ' ) 

KK*-11 
260  KK*KK-»-12 

L=MIN0(KK+11,N) 

DO  270  1*1, N 

270  WRITE(12,110)  (K(I, J) , J=KK,L) 

IF(L.LT.N)GO  TO  260 
C  DO  60  1*1,12 

C  WRITEd2,101)  (K(I,J),J=1,9) 

C  60  CONTINUE 

C  WRITE (12, 102) 

C  DO  61  1*1,12 

C  WRITE(12,101)  (K(I,J) ,J*10.12) 

C  61  CONTINUE 

C  WR1TE(12,102) 

C  DO  62  1*1,21 

C  WRITE(12,104)  (K(I,J),J*19,21) 

C  62  CONTINUE 

C.  •  . 

C...  FORMATS 

C.  .  . 

101  FORMAT (9F12. 8) 

C102  FORMATdH  ) 

102  FORMAT (IX) 

103  FORMAT (A5) 

104  FORMAT (3F12. 8) 

105  FORMAT(/) 

110  FORMAT (12E15. 6) 
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c 

c 

C  40 

c 

c 

c 

C  41 
c. . . 
c. . . 
c. . . 


50 

c 

c 

c 

C51 

C 

C 

C 

C52 


111  FORMAT (12F10.6) 

112  FORMAT (9F10. 6) 

114  FORMAT (12E15. 6) 

121  FORMAT ( IX, 3E20. 12) 

RETURN 

END 

SUBROUTINE  MATCALC(A, B,N.M) 

C... 

C. . .  THIS  SUBROUTINE  WILL  DETERMINE 

C...  (1)  DET  OF  A 

C...  (2)  INVERSE  OF  A 

C...  (3)  SOLVE  A  SYSTQ<  OF  E(2UATIC»1S 

C. . .  BASED  ON  THE  VALUE  OF  THE  PARAMETER  INDEX 

C...  IF  INDEX  EQUALS  (0.1,-1)  THE  OPTION  SHOWN  ABOVE  WILL  BE  DETERMINED 
C...  A(N,M)  a  TH'^  AUOfENTED  MATRIX 

C...  B(N,N}  a  ORIGINALLY  THE  NxN  IDENTITY. .  .THE  INVERSE  MATRIX  FINALLY 
C . .  .  THE  METHOD  USED  IS  GAUSSIAN  ELIMINATION  WITH  PIVOTING 
C.  .  . 

IMPLICIT  REAL* 8  (A-H,0-2) 

DIMENSION  A(N,M)  ,B(N,M) 

CCAfMON  NSYS,  INDEX.  DET 

SIGNsl 

MARK  «0 

NMIsN-l 

NN=2*N 

NPLSYaN-i-NSYS 

IF(INDEX.LE.O)  GO  TO  2 

DO  1  1=1, N 

DO  1  J=1,N 

1  A(I,N-fJ)=B(I,  J) 

NPLSY=NN 

2  CONTINUE 

DO  10  1=1, NMI 

C. .  . 

C . . .  FROM  HERE  TO  STATEMENT  4  THE  PROGRAM  PICKS  UP  THE  PIVOT 
C.  .  . 

MAXal 

AMAX=ABS(A(I,I)) 

K=I 

3  K=K-»-l 

IF(ABS(A(K,I) ) .LE.AMAX)  GO  TO  4 
MAX=K 

AMAXaABS(A(K,I) ) 

4  IF(K.NE.N)  GO  TO  3 
IF(MAX.EQ.I)  GO  TO  6 

C.  . . 

C. . .  THE  NEXT  SEQUENCE  INTERCHANGES  ROWS 
C.  .  . 

L=I-1 
J  L=L-i-l 

TEMP=A(I,L) 

A(I,L)=A(MAX,L) 

A(MAX,L}=TEMP 
IF(L.LT.NPLSY)  GO  TO  5 
SIGN=-SIGN 

6  J»I 

7  J=J+1 

IP(A(J,I) .EQ.0.0)  GO  TO  9 
CONSTa-A(J,I)/A(I,I) 
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8 


L-I-1 

L-lrfl 

A(J,  L)  «A(J,  L) -t-ACI.  L)  *CONST 
IF(L.NE.NPLSY)  GO  TO  8 

9  C(»n'lNUE 
IF(J.NE.N)  GO  TO  7 

10  CONTINUE 
TEMPsl 

DO  11  1«1.N 

IF(A(I.Z) .EQ.0.0)  GO  TO  12 

11  TEHP>TEMP*A(I,Z) 

OET>SI(ai*TQfP 
GO  TO  13 

12  HARKsl 
DET«0.0 

13  IF(ZNDBX.EQ.O)  GO  TO  21 
IF(MARK.NE.l)  GO  TO  15 
NRZTE(6,14) 

C. . . 

C...  FORMATS 
C.  . . 

14  F0RMAT(///2X,21HMATRIX  A  IS  SINGULAR.) 

GO  TO  21 

15  NlsN-t-1 
C. . . 

C...  HERE  THE  PROGRAM  CARRIES  OUT  BACK  SUBSTITUTION 
C. .  . 

DO  20  laNl.NPLSY 
KaN 

16  B(K,I)aA(K,I) 

IF(K.EQ.N)  GO  TO  18 
JaK 

17  JaJ+1 

B(K. I)aB{K. I)-A(K, J) *B(J, I) 

IF(J.NE.N)  GO  TO  17 

18  B(K,I)aB(K,I)/A(K,K) 

IF(K.EQ.l)  GO  TO  19 
KaK-1 

GO  TO  16 

19  CONTINUE 

DO  20  Lal.N 

20  A(L,I)aB(L,I) 

21  RETURN 
END 

SUBROUTINE  OUTPAK  (MATRIX,  NROW,MMATR,NCTL,NOUT) 

C . VERSION  a  09/05/73/03 

C . 

c 

C  OUTPAK  PRINTS  A  REAL*8  SYMMETRIC  MATRIX  STORED  IN  ROW-PACKED  LOWER 
C 

C  TRIANGULAR  FORM  (SEE  DIAGRAM  BELOW)  IN  FORMATTED  FORM  WITH  NUMBERED 
C 

C  ROWS  AND  COLUMNS.  THE  INPUT  IS  AS  FOLLOWS; 


C 

C  MATRIX(*) . PACKED  MATRIX 

C 

C  NROW . NUMBER  OF  ROWS  TO  BE  OUTPUT 

C 

C  NCTL . CARRIAGE  CONTROL  FLAG:  1  FOR  SINGLE  SPACE. 

-21- 


C  2  FOR  DOUBLE  SPACE, 

C  3  FOR  TRIPLE  SPACE. 

C 

C  NOOT . UNIT  NUMBER  FOR  OUTPUT 

C 


C  THE  MATRIX  ELEMENTS  ARE  ARRANGED  IN  STORAGE  AS  FOLLOWS: 
C 

C  1 

C  2  3 


c 

4 

5 

< 

c 

7 

8 

S 

10 

c 

11 

12 

13 

14 

15 

c 

16 

17 

18 

19 

20 

21 

c 

22 

23 

24 

25 

26 

27 

C 

C  AND  SO  ON. 

C 

C  OUTPAK  IS  SET  UP  TO  HANDLE  6  COLUMNS/PAGE  WITH  A  SP20.14  FORMAT 
C 

C  FOR  THE  COLUMNS.  IF  A  DIFFERENT  NUMBER  OF  COLUMNS  IS  REQUIRED,  CHANGE 
C 

C  FORMATS  1000  AND  2000,  AND  INITIALIZE  KCOL  WITH  THE  NEW  NUMBER  OF 
C 

C  COLUMNS. 

C 

C  AUTHOR:  NELSON  H.F.  BEEBE,  QUANTUM  THEORY  PROJECT,  UNIVERSITY  OF 
C  FLORIDA,  GAINESVILLE 

C 

C . 

REAL*8  MATRIX, COLUMN 
INTEGER  BEGIN. ASA, BLANK, CTL 
DIMENSION  MATRIX ( NMATR ) , ASA ( 3 ) 

DATA  KCOL/5/,  COLUMN/ 8HCOLUMN  /,  ASA/4H  ,  4H0  ,  4H-  /, 

&  BLANK/4H  /,  ZERO/O.D+00/ 

CTL  s  BLANK 

IF  ( (NCTL.LE.3) .AND. (NCTL.GT.O))  CTL  *  ASA{NCTL) 

C 

C  LAST  IS  THE  LAST  COLUMN  NUMBER  IN  THE  ROW  CURRENTLY  BEING  PRINTED 
C 

LAST  s  MINO(NROW,KCOL) 

C 

C  BEGIN  IS  THE  FIRST  COLUMN  NUMBER  IN  THE  ROW  CURRENTLY  BEING  PRINTED. 

C 

NCOLsl 

C . BEGIN  NON  STANDARD  DO  LOOP. 

BEGINsl 
1050  NCOL  «  1 

WRITE  (NOUT,1000)  (1,1  s  BEGIN,LAST) 

DO  40  K  s  BEGIN, NROW 
KTOTAL  »  (K*(K-l))/2  ♦  BEGIN  -  1 
C  DO  10  I  «  l.NCOL 

C  GO  TO  20 

C  IF  (MATRIX (KTOTAL+ I)  .NE.  ZERO)  GO  TO  20 
C  10  CONTINUE 
C  GO  TO  30 

20  WRITE  (NOUT,2000)  CTL, K, (MATRIX(I+KTOTAL) , I»l, NCOL) 

30  IF  (K  .LT.  (BEGIN-i-KCOL-l))  NCOL  s  NCOL  *  1 
40  CONTINUE 

LAST  «  MINO(LAST-i-KCOL,NROW) 

BEGINsBEGIN-fNCOL 
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IP  ( BEGIN. LE-NROW)  GO  TO  lOSO 
1000  FORMAT  (/12X, 4 (IIX, 14. 5X) , (IIX. 14) ) 
2000  FORMAT  (Al, 4X, 14. 2X. 5D20 .12) 

RETURN 

END 
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IL2.3  Program  inatiniilt2J 

The  third  program  is  matmult2i.  This  program  also  carries  out  die  transfonnatitm 

F=B^KB-*  (8) 

and  also  determines  K  from  F.  The  program  allows  input  d  the  scaling  factors,  die  Qj’s,  to  scale 
die  force  constant  matrix  in  internal  coordinatBS,  F,  and  converts  the  scaled  F  to  a  scaled  force 
constant  matrix  in  Cartesian  coordinates.  K.  This  scaled  K  is  used  as  iiqiut  to  die  CADPAC 
program  to  carry  out  a  VCD  calculation  of  allowed  frequencies  of  vibration  and  rotatirmal 
strengths.  In  addition .  the  parameter  NAT  must  be  changed  for  each  molecule  considered. 

On  the  following  pages  a  listing  of  the  FORTRAN  program  nuUmult2.f  is  ^en. 
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iDatioult2.£ 


C. . . 
C. . . 
C. . . 


C. . 
C. . 
C. . 
C. . 
C. . 
C. . 

c. . 
c. . 
c, . 
c. . 
c. . 
c. . 
c. . 
c. . 
c. . 
c. , 
c. . 
c. . 
c. . 
c. , 
c. . 

c. . 
c. . 
c. . 


10 

c. .  . 
c. .  . 
c. .  . 


501 


502 

C.  .  . 


PROGRAM  FCMATRIX 

PUNCHES  F.C.M  TO  FORTRAN  UNIT  7 

PARAMETER  (NATsl6 .  MM«3  ‘NAT-fi .  N«3  *MAT.  MM1«2  *MM,  NAT3  «N,  NDIM«MAT3 ) 
IMPLICIT  REALMS  (A-H,0>Z) 

REAL*8  K(N.N).M(N.N).KMEW(M,N) 

REAL  MC.MO.MH 

DIMENSION  B(MM,N).BM(MM,N)  ,BP(N.MM),P(N,M),TEST(MM,MM) , 

1  BMBP(MM,MH)  ,BMBPI(»0<.M<).BPLM(MM,N}.PROD(MM,N},BPLM1(N,M(}. 

2  TESTA  (MM,  MM) .  AA(MM.MMM}  ,BB(MM,)«M)  .FNEW(MM,MM) , 

3  TITLE  ( 10 ) .  C  ( 3 ,  NAT) .  GRAD  ( 3 .  NAT) ,  FCM(NDIM,  NDDf) ,  PROD2  (MM,  N) , 

4  Q(15) 

CCB040N  NSYS,  INDEX,  DET 

READ(8,105)  (TITLE(I),I*1,9) 

WRITE(7,105)  (TITLE(I) ,Ial,9) 

READ(8,10€) 

WRITE(7,106) 

READ(8,107)  (C(l,I),C(2,I),C(3,I),Isl,NAT) 

WRITE(7,107)  (C(1,I),C(2,I),C(3,1),I«1,NAT) 

R£AD(8,10d) 

WRITE (7, 108) 

READ(8,107)  (GRAD(1,I) ,GRAD(2, I) ,GRAD(3, 1) , I*1,NAT) 

WRITE(7.107)  (GRAD(1,I),GRAD{2,I),GRAD(3,I),I«1,NAT) 

READ(8,109) 

WRITE (7, 109) 

GET  MATRICES  B  AND  K 

CALL  BKMATR(MM,N,B,K) 

ADJOINT  OF  B 

DO  10  I«1,N 
DO  10  Jsl,MM 

BP(I,J).B(J,I) 

CONTINUE 

DETERMINE  PRODUCT  OF  B  M  BP  MATRICES 

MC«12.01 
MO«16.00 
MH«1.008 
DO  501  I>1,N 
DO  501  Jsl,N 
M(I, J)*0, 

CONTINUE 
DO  502  I«1,N 
M(I,I)«1. 

CONTINUE 
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DO  11  I«1,MM 
DO  11  J«1,N 

SUM«0. 

DO  11  L«1,N 

SUMsSUM-fB  ( I .  L)  *M (L,  J) 

BM(Z,  J)*SX]M 

11  CONTINUE 
C. . . 

C. . .  DETERMINE  PRODUCT  OF  B  M  BP  MATRICES 

DO  511  I«1.MM 
DO  511  J>1,MM 
SUM>0. 

DO  511  L«1,N 

SUMsSUM-t-BMd.L)  *BP(L,  J) 
BMBP<I, J}«SUM 

511  CONTINUE 

OPEN ( 3 , FILES ' TESTADZ . OUT ' ) 

OPEN ( 23 , FILES ' BMBP .mat ' } 

WRITE{3,*)  'BMBP' 

DO  540  Isl.MM 

WRITE(3,115)  (BMBP(I,J).Js1,MM) 
WRITE(23,115)  (BMBP(I,J),Js1,MM) 

540  CONTINUE 

C.  .  . 

C . . .  DETERMINE  INVERSE  OF  B  H  BP 

DO  12  Isl.MM 
DO  12  Jsl.MM 

AA(I. J)sBMBP(I,J) 

12  CONTINUE 
NSYSsO 
INDEXsl 

DO  221  Isl.MM 
DO  221  JsMM^l.MMM 
AA(I,J)sO. 

IF( (J-MM) .EQ.I)  AA(I.J)=1. 

221  CONTINUE 

DO  222  Isl.MM 
DO  222  Jsl.MM 
BB(I. J)s0. 

IF(I.EQ.J)  BB(I.J)s1. 

222  CONTINUE 

CALL  MATCALC(AA.BB,MM.MMM) 

WRITE(6,*)  'DETA  s'.DET 

C.  .  . 

C...  SET  BMBPI  MATRIX 

C. . . 

DO  191  Isl.MM 
DO  191  Jsl.MM 

BMBPI ( I . J) sAA ( I , J+MM) 

191  CONTINUE 
WRITE(3.*)  'BMBPI' 

DO  192  Isl.MM 

WRITE(3.116)  (BMBPI(I.J) .Jsl.MM) 

192  CONTINUE 
C. . . 

C...  DETERMINE  TESTA  MATRIX 
C.  .  . 


DO  302  Isl.MM 
DO  302  Jsl.MM 
SUMsO. 
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DO  302  L«1,1M 

SUM*SUM-i-BMBPI  (X.L)  *BMBP(L,  J) 

TESTAd,  J)>SUM 
302  CONTINUE 

WRITE(3,*)  'TESTA' 

DO  340 

WRITE(3,111)  (TESTAd, 

340  CONTINUE 

C  NRITE(3,102) 

C  DO  341  I«l,15 

C  WRITE(3,112)  (TESTAd,J),J«13,15) 

C  341  CONTINUE 
C. . . 

C.  . .  DETERMINE  BPLM  MATRIX 

C.  . . 

DO  13  Z-l.MM 
DO  13  J«1,N 

SUMsO. 

DO  13  L«1,MM 

SUMsSUM-i-BMBPI  (I ,  L)  *BM  (L.  J) 

BPLMd,  J)«SUM 

13  CONTINUE 
C.  .  . 

C.  .  .  DETERMINE  TEST  MATRIX 

C.  .  . 

DO  202  Isl.MM 

DO  202  Jsl.MM 

SUMsO. 

DO  202  L=1,N 

SUMs:SUM4-BPLM  ( I .  L)  *BP  (L,  J) 

TESTd,  J)*SUM 

202  CONTINUE 

OPEN ( 2 , FILE= ' TESTDZ . OUT ' ) 

DO  240  Isl,MM 

WRITE{2,111)  (TESTd,J)  .J*1,MM) 

240  CONTINUE 

C  WRITE(2,102) 

C  DO  241  1=1,15 

C  WRITE(2,112)  (TESTd,  J)  ,J=13,15) 

C  241  CONTINUE 

C.  .  . 

C . . .  DETERMINE  TRANSPOSE  OF  BPLM  MATRIX 
C.  .  . 

DO  14  1=1, N 
DO  14  J«1,MM 

BPLMI (I , J ) «BPU<  ( J , I ) 

14  CONTINUE 
C. . . 

C...  DETERMINE  PRODUCT  OF  BPUl  K  MATRICES 
C.  .  . 

DO  20  1=1, MM 
DO  20  J«1,N 

SUM=0. 

DO  20  L«1,N 

SUM=SUM+BPLMd,L)  *K(L,  J) 

PRODd,J)*SUM 

20  CONTINUE 

C.  .  . 

C . .  .  DETERMINE  PRODUCT  OP  BPLM  K  BPLMI  MATRICES  TO  GIVE 

C...  Pd,J)  MATRIX 
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o  rt  n 


C. . . 

DO  30 

DO  30  Jsl.MM 
SUM«0. 

DO  30  L«1,N 

SUMeSUM-^PROD ( I .  L)  *BPLMI  (L,  J) 

F(I,J)>SUM 

C. . 

C. . .  IN  ORDER  TO  CONVERT  TO  UNITS  OP  MDYNE/A 
C. . .  INSERT  THE  FOLLOWING  STATEMENT 

C.  .  . 

C  F(I.J)«15.57*F(I.J) 

30  CONTINUE 

C. . . 

C...  F(I,J)  MATRIX 

C. . . 

OPEN (1, FILE* 'FORCE. OUT' ) 

WRITE (1.*)  '  FORCE  CONSTANT  MATRIX' 

WRITE (1.*)  '  (INTERNAL  COORDINATES  -  UNITS  OF  HARTREES/A) ' 
WRITEd.*)  ' 

DO  40  Isl.MM 

WRITEd,  111)  (Fd,  J) ,  J=1,MM) 

40  CONTINUE 

C  WRITEd,  102) 

C  DO  41  1=1,15 

C  WRITEd, 112)  (F(I,  J) ,  J=13,15) 

C  41  CONTINUE 

SCALING  FACTORS  Q(I)  INPUT  HERE 

Qd)=  0.958 
Q(2)=  0.958 
Q(3)=  0.958 
Q(4)=  0.907 
Q(5)=  0.772 
Q(6)=  0.863 
Q(7)=  0.931 
Q(8)=  0.845 
Q{9)*  0.907 
Qd0)=  0.863 
Qdl)=  0.845 
Qd2)=  0.907 
Qd3)=  0.863 
Qd4)=  0.863 
Qd5)=  0.845 
Qd6)=  0.923 
0(17)*  0.923 
Q(18)«  0.914 
0(19)*  0.904 
0(20)*  0.901 
0(21)*  0.903 
0(22)*  0.946 
0(23)*  0.902 
0(24)*  0.901 
0(25)*  0.946 
0(26)*  0.902 
0(27)*  0.902 
0(28}*  0.901 
0(29)*  0.946 
0(30)*  0.900 


601 

1230 

1231 


807 


701 
C. . . 
C.  . . 
C. .  . 


C 

602 


603 

C. .  . 
C 


660 
C. . . 


Q(31)«  0.932 
Q(32)«  0.916 
Q(33)>  0.900 
Q(34)«  0.900 
g(35)«  1.093 
Q(36)«  0.921 
Q(37)s  0.910 
g(38)«  1.093 
Q(39)«  0.921 
Q(40)m  0.921 
Q(41)«  0.910 
Q(42)>  0.984 

NEW  P  MATOIX,  FNEW 

OPEN ( 5 1 . FILE- ' FNEW . OUT ' ) 

WRITE (51,*)  'FNEW  IN  UNITS  OP  HARTREES/BOHR' 

DO  601  1-1. MM 
DO  601  J«1,MM 

FNEW(I,J)«SQRT(Q(I)*Q(J))*F{I,J) 

IPd.EQ.J)  FNEW(I,J)«Q(I)*F(I,J) 

CONTINUE 

FORMAT(15E15.6) 

FORMAT(I5.E15.6) 

DO  807  1-1, MM 

WRITE (51, 1230)  (FNEW(I, J) , J-1, I) 

CONTINUE 
WRITE (51, 102) 

WRITE (51, 102) 

WRITE(51,103) 

WRITE (51, 102) 

DO  701  I-l.MM 
WRITE (51, 1231)  I,Q(1) 

CONTINUE 

NEW  K  MATRIX,  KNEW 

DO  602  I-1,MM 
DO  602  Jsl,N 
SUM-0. 

DO  602  L-1,MM 

SUM«SUM-*-FNEW(I,L)  *B(L,  J) 

SUM-SUM-t-FNEW  ( I ,  L)  / 15 . 57  *B  (L,  J) 

PROD2(I,J)«SUM 
CONTINUE 
DO  603  I-1,N 
DO  603  J«1,N 
SUM-0. 

DO  603  L«1,MM 

SUM-SUM-t-BP  ( I ,  L)  *PR002  (L,  J) 

KNEW(I,J)«SUM 
PCM(I, J)*KNEW(I,  J) 

CONTINUE 
DO  660  I-1,NAT 

READ(8,107)  (FCM( J, 1*3-2), FCM(J, 1*3-1), FCM(J,  1*3 ) ,J-1,NAT3) 
OPEN { 7 1 , FILE- ' KMATNEW . OUT ' ) 

WRITE(71,107)  (FCM(J,1*3-2),FCM(J,I*3-1),FCM(J, 1*3) ,J*1,NAT3) 
CONTINUE 
FORMATS 
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c. .  . 

105  FORMAT (IX, 9 AS) 

106  FORMATdX,  'GEC»«ETRY' ) 

107  FORMATdX,  3E20. 12) 

108  FORMATdX,  'GRADIENT' ) 

109  FORMATdX, 'CARTESIAN  SECOND  DERIVATIVES  (UNPROJECTED)') 

102  FORMATdX) 

103  FORMAT(4X, 'I',11X, 'Q(I) ') 

111  FORMAT (12F12. 6) 

112  FORMAT(9F12.6) 

115  FORMAT (15F12. 6) 

116  FORMAT (15E15. 6) 

STOP 

END 

SUBROUTINE  BKMATR(M,N,  B,K) 

IMPLICIT  REALMS  (A-H,0-Z) 

REAL*8  K 

DIMENSION  B(M,N) ,K(N,N) 

NAT»N/3 

MMsN-6 

WRITE(*,*)  'NAT  =  ',NAT 
WRITE(*,*)  'MM  *  ',MM 
C , . .  B  MATRIX 

OPEN(21,FILE='BMAT.IN' ) 

KK*-11 

160  KK=KK+12 
L=MIN0(KK+11,N) 

DO  170  1=1, MM 

170  READ(21,114)  (8(1, J) , J=KK,L) 

IF(L.LT.N)GO  TO  160 

C  DO  30  1=1, MM 

C  READ(21,114)  (B(I,J),J=1,12) 

C  30  CONTINUE 

C  DO  31  1=1,15 

C  READ(21,114)  (B(I,J),J=13,21) 

C  31  CONTINUE 

OPEN(22,FILE='BMAT.OUT' ) 

KK=-11 

161  KK=KK+12 
L=MIN0(KK+11,N) 

DO  171  1=1, MM 

171  WRITE(22,114)  (B(I, J) , J=KK,L) 

IF(L.LT.N)GO  TO  161 

C  DO  40  1=1, MM 

C  WRITE(22,111)  (B(I,J),J=1,12) 

C  40  CONTINUE 

C  WRITE(22,102) 

C  DO  41  1=1,15 

C  WRITE(22,112)  (B(I, J) , J=13,21) 

C  41  CONTINUE 

C.  .  . 

C . . .  K  MATRIX 

C.  .  . 

OPEN (1 1 , FILE= ' KMAT . IN ' ) 

DO  50  1=1, NAT 

READdl,121)  (K(  J,  1*3-2),  K(J,  1*3-1), K(J,  1*3),  J=1,N) 

50  CONTINUE 

C  READ(11,105) 

C  DO  51  1=1,21 


-30- 


C  R£AD(11,101)  (K(I, J),J-10,18) 

C51  CONTINUE 

C  It£AO(ll,10S) 

C  DO  52  Isl,21 

C  READ(11,104)  (Kd.  J)  ,J«19.21) 

C52  CONTINUE 

OPEN ( 12 . PILE* ' KMXT .OUT' ) 

DO  60  Iml.S 

WRITE(12,101)  (K(I.J).J.1.9) 

60  CONTINUE 
WRITE (12. 102) 

DO  61  I«1,N 

WRITE(12,101)  (K(I.J),J«10,12) 

61  CONTINUE 

C  WRITE (12, 102) 

C  DO  62  l3l,21 

C  WRITE(12,104)  (K(I, J) , J-19,21) 

C  62  CONTINUE 

C.  .  . 

C . . .  FORMATS 

C.  .  . 

101  FORMAT (9F12. 8) 

C102  FORMAT (IH  ) 

102  FORMAT (IX) 

103  FORMAT (AS) 

104  FORMAT(3F12.8) 

105  FORMAT(/) 

111  FORMAT (12F10. 6) 

112  FORMAT (9F10. 6) 

114  FORMAT (12E15. 6) 

121  FORMAT { IX, 3E20. 12) 

RETURN 

END 

SUBROUTINE  MATCALC (A, B,N,M) 

C.  .  . 

C . . .  THIS  SUBROUTINE  WILL  DETERMINE 

C. . .  (1)  DET  OF  A 

C.  .  .  (2)  INVERSE  OF  A 

C...  (3)  SOLVE  A  SYSTEM  OF  EQUATIONS 

C .  .  .  BASED  ON  THE  VALUE  OF  THE  PARAMETER  INDEX 

C...  IP  INDEX  EQUALS  (0,1,-1)  THE  OPTION  SHOWN  ABOVE  WILL  BE  DETERMINED 
C...  A(M,M)  s  THE  AUOIENTED  MATRIX 

C...  B(N,N)  =  ORIGINALLY  THE  NxN  IDENTITY .. .THE  INVERSE  MATRIX  FINALLY 
C...  THE  METHOD  USED  IS  GAUSSIAN  ELIMINATION  WITH  PIVOTING 
C. . . 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  A(N,M) ,B(N,M) 

COMMON  NSYS, INDEX, DET 

SIGN=1 

MARK  sO 

NMIsN-1 

NN=2*N 

NPLSYaN+NSYS 

IF(INDEX.LE.O)  GO  TO  2 

DO  1  1*1, N 

DO  1  J=1,N 

1  A(I,N+J)*B(I, J) 

NPLSYaNN 

2  CONTINUE 


DO  10 

C. .  . 

C.  . .  FROM  HERE  TO  STATQfENT  4  THE  PROGRAM  PICKS  UP  THE  PIVOT 
C.  .  . 

MAXsI 

AMAX-ABS(A(I,I) ) 

K«I 

3  K>.K4-1 

IF(ABS(A(X,I)} .LE.AMAX)  GO  TO  4 
MAX«K 

AMAXsABS(A(K,I)) 

4  IF(K.NE.N)  GO  TO  3 

IF(MAX.EQ.I)  GO  TO  6 

C.  . . 

C. . .  THE  NEXT  SEQUENCE  INTERCHANGES  ROWS 
C.  . . 

L»I-1 

5  LsL-t- 1 
TEMP=AII.Li 
A(I.li}«A(MAX.L) 

A(MAX,L)«TEMP 
IF(L.LT.NPLSY)  GO  TO  5 
SIGNs-SIGN 

6  J=1 

7  JsJ+1 

IF(A(J,I) .EQ.0.0)  GO  TO  9 
CONST=-A(J,I)/A(I,I) 

L=l-1 

8  LsL-fl 

A(J,L)»A(J,L)+A(I,L)*CONST 
IF(L.NE.NPLSY)  GO  TO  8 

9  CONTINUE 
IF(J.NE.N)  GO  TO  7 

10  CONTINUE 
TEMPsl 

DO  11  1=1, N 

IF(A(I,I) .EQ.0.0)  GO  TO  12 

11  TEMP=TEMP*A(I, I) 

DET=SIGN*TEMP 
GO  TO  13 

12  MARK=1 
DET=0 . 0 

13  IF(INDEX.EQ.O)  GO  TO  21 

IF(MARK.NE.l)  GO  TO  15 
WRITE(6,14) 

C.  .  . 

C...  FORMATS 

C.  . . 

14  FORMAT (///2X,21HMATRIX  A  IS  SINGULAR.) 

GO  TO  21 

15  N1=N+1 
C. .  . 

C . . .  HERE  THE  PROGRAM  CARRIES  OUT  BACK  SUBSTITUTION 

C.  .  . 

DO  20  I=N1,NPLSY 
K=N 

16  B(K,I)=A(K.I) 

IF(K.EQ.N)  GO  TO  18 
J*K 

17  J=5+l 
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B(K.  I)-A(K,  J)  *B(J.  I) 

IF(J.NE.N)  GO  TO  17 
B(K,I)>B(K.I)/A(K,K) 
IF(K.EQ.l)  GO  TO  19 
K«K-1 
GO  TO  16 

19  CONTINUE 

DO  20  L«1,N 

20  A(L,I)«B(L,I) 

21  RETURN 
END 
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n.2.4  Program  sinnplex.f 

Another  program  named  simplex.f  was  generated  from  one  of  the  programs  of  the  programs 
of  McIntosh  and  Peterson  [1].  This  program  allows  the  scaling  factors  to  be  determined  by  a 
best  fit  to  a  set  of  inputted  experimental  frequencies.  This  approach  was  investigared;  however, 
the  approach  used  in  this  st^y  was  to  determine  a  Qj  1^  comparison  of  the  diagonal  force 
constants  at  the  6*3  IG*  level  F-CHF)  and  F-(M^).  The  simplex  approach  may  be  the  better 
approach  and  should  be  given  serious  consiaamion  for  scaling  pnx^utes  to  be  studied  in  the 
future. 
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Table  1.  Data  file  for  R-glyceraldehyde  which  is  used  as  input  to  die  FORTRAN  program 
bmatf.  The  program  bmatf  determines  the  transfcmnation  matrix  B  defined  1^  R  s  B  q  where 
R  is  a  column  vector  of  the  Cartesian  coordinates. 


C3H603  -  glycarald^hyd*  [hf/C-31g*] 
12  0 


-2.026663 

- 

■1.245592 

0.698964 

C-1 

-0.075542 

0.791716 

1.181527 

C-2 

2.590801 

• 

-0.265346 

0.987662 

C-3 

-1.840047 

- 

-2.980089 

1.800025 

H-4 

-3.709799 

- 

-0.957621 

- 

-0.760258 

0-5 

-0.389204 

1.489681 

3.097911 

H-6 

-0.333856 

2.731754 

- 

-0.577374 

0-7 

-1.950507 

2.587079 

- 

-1.351823 

H-8 

3.928350 

1.233993 

1.405427 

H-9 

2.869625 

- 

-1.763088 

2.355118 

H-10 

3.009882 

- 

-1.303036 

- 

-1.399301 

0-11 

2.720395 

-0.021029 

• 

-2.620120 

H-12 

1 

1 

2 

0 

0 

0 

0 

1-2  bond  stretch 

1 

2 

3 

0 

0 

0 

0 

2-3  bond  stretch 

1 

1 

4 

0 

0 

0 

0 

1-4  bond  stretch 

1 

1 

5 

0 

0 

0 

0 

1-5  bond  stretch 

1 

2 

6 

0 

0 

0 

0 

2-6  bond  stretch 

1 

2 

7 

0 

0 

0 

0 

2-7  bond  stretch 

1 

7 

8 

0 

0 

0 

0 

7-8  bond  stretch 

1 

3 

9 

0 

0 

0 

0 

3-9  bond  stretch 

1 

3 

10 

0 

0 

0 

0 

3-10  bond  stretch 

1 

3 

11 

0 

0 

0 

0 

3-11  bond  stretch 

1 

11 

12 

0 

0 

0 

0 

11-12  bond  stretch 

2 

3 

2 

1 

0 

0 

0 

bond  angle  bend  3-2-1 

2 

4 

1 

2 

0 

0 

0 

bond  euigle  bend  4-1-2 

2 

5 

1 

2 

0 

0 

0 

bond  angle  bend  5-1-2 

2 

6 

2 

1 

0 

0 

0 

bond  angle  bend  6-2-1 

2 

7 

2 

1 

0 

0 

0 

bond  angle  bend  7-2-1 

2 

8 

7 

2 

0 

0 

0 

bond  angle  bend  8-7-2 

2 

9 

3 

2 

0 

0 

0 

bond  angle  bend  9-3-2 

2 

10 

3 

2 

0 

0 

0 

bond  angle  bend  10-3-2 

2 

11 

3 

2 

0 

0 

0 

bond  angle  bend  11-3-2 

2 

12 

11 

3 

0 

0 

0 

bond  angle  bend  12-11-3 
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Table  1.  Data  file  for  R-glyceraldehyde  wMch  is  used  as  input  to  the  FORTRAN  program 
Imatf.  The  program  bmatf  determines  the  transformation  nutiix  B  defined  by  R  «  B  q  v^iere 
R  is  a  column  vector  of  the  Cartesian  coordinates.  (CONTINUED) 


dihadral  «ngl*  4>l-2-3 
dihadral  angla  5-1-2-3 
dihedral  angle  6-2-1-4 
dihedral  angle  7-2-1-4 
dihedral  angle  8-7-2-1 

dihedral  angle  9-3-2-1 
dihedral  angle  10-3-2-1 
dihedral  tmgle  11-3-2-1 
dihedral  angle  12-11-3-2 
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